% fig6a.m
clear; close all; clc;

%% 1. 模型参数
alpha = 5;
sigma = 0.8;
mu    = 0.1;
epsilon = 0.1;

N_total = 12000;  % 总迭代次数
N_cut   = 100;  % 去除瞬态

% 扫描范围: phi0 in [-4, 6], k in [0, 0.4]
phi_min = -4; phi_max = 6; Nphi = 301;
k_min   =  0; k_max   = 0.4; Nk   = 301;

phi_values = linspace(phi_min, phi_max, Nphi);
k_values   = linspace(k_min,   k_max,   Nk);

%% 2. 预分配分类矩阵 cat_map
cat_map = zeros(Nphi, Nk);

%% 3. 双参数扫描
parfor i_phi = 1:Nphi
    phi0_now = phi_values(i_phi);
    
    tmp_row = zeros(1, Nk);
    for i_k = 1:Nk
        k_now = k_values(i_k);
        
        % 调用分类函数
        catLabel = classifyDynamics(phi0_now, k_now, alpha, sigma, mu, epsilon, N_total, N_cut);
        
        tmp_row(i_k) = catLabel;
    end
    
    cat_map(:, i_phi) = tmp_row;
end

%% 4. 绘制分类结果
figure;
imagesc(k_values, phi_values, cat_map);
axis xy;
xlabel('k');
ylabel('\phi_0');
title('(a) 动力学分类');


cmap = [
    1.0, 0.0, 0.0;  % 红色: HC (Hyperchaos, 超混沌)
    1.0, 0.5, 0.0;  % 橙色: CH (Chaos, 混沌)
    1.0, 1.0, 0.0;  % 黄色: QP (Quasi-Periodic, 准周期)
    1.0, 0.5, 1.0;  % 粉色: MP (Multiple Periods, 多周期)
    0.8, 0.5, 0.0;  % 棕橙色: P7 (周期7)
    0.5, 0.4, 0.0;  % 棕色: P6 (周期6)
    0.0, 1.0, 1.0;  % 青色: P5 (周期5)
    0.0, 0.7, 0.7;  % 深青色: P4 (周期4)
    0.0, 0.6, 0.2;  % 橄榄绿色: P3 (周期3)
    0.0, 0.3, 0.5;  % 蓝绿色: P2 (周期2)
    0.0, 0.0, 1.0;  % 蓝色: P1 (周期1)
];
colormap(cmap);

% colorbar 范围固定在 [1, 11]
clim([1 12]);
cb = colorbar;
cb.Ticks = 1.5:1:11.5;  % 中间值
cb.TickLabels = {'HC','CH','QP','MP','P7','P6','P5','P4','P3','P2','P1'};
title('(a) 动力学分类');


function catLabel = classifyDynamics(phi0, k, alpha, sigma, mu, epsilon, N_total, N_cut)
% classifyDynamics 改进后的分类函数，结合前两个 Lyapunov 指数 (由 LEs.m 计算)
% 和周期判定对系统状态进行分类。输出分类标签：
%        1 -> HC (Hyperchaos, 超混沌)
%        2 -> CH (Chaos, 混沌)
%        3 -> QP (Quasi-Periodic, 准周期)
%        4 -> MP (Multiple Periods, 多周期)
%        5 -> P7 (周期7)
%        6 -> P6 (周期6)
%        7 -> P5 (周期5)
%        8 -> P4 (周期4)
%        9 -> P3 (周期3)
%       10 -> P2 (周期2)
%       11 -> P1 (周期1)

    % 1. 运行 Rulkov 模型获得时间序列及归一化同步误差 E_sim
    [X1_save, ~, E_sim] = Rulkov(alpha, sigma, mu, k, phi0, N_total, epsilon, N_cut);
    idx = (N_cut+1):N_total;
    X1_eff = X1_save(idx);
    
    % 2. 先根据同步误差初步分类（同步误差低时认为是强同步周期 P1；中等同步误差判为混沌 CH）
    if E_sim < 0.05
         catLabel = 11; % P1: 强同步
         return;
    elseif E_sim < 0.15
         catLabel = 2;  % CH: 混沌
         return;
    end
    
    % 3. 调用 LEs.m 计算前两个 Lyapunov 指数
    [LE1, LE2] = LEs(alpha, sigma, mu, k, phi0, epsilon);
    
    % 4. 利用 Lyapunov 指数判断超混沌和准周期
    %    此处阈值为示例，建议根据仿真结果调整
    if (LE1 > 0.2) && (LE2 > 0.1)
         % 两个正指数均较大，表明系统为超混沌
         catLabel = 1; % HC: 超混沌
         return;
    elseif (LE1 > 0.05) && (LE1 <= 0.2)
         % 如果 LE1 较小但正，判为混沌 CH（前面 E_sim 已排除极低同步状态）
         catLabel = 2; % CH: 混沌
         return;
    end
    
    % 5. 对剩余情况，利用周期检测进一步分类
    [p, cv] = Period(X1_eff);   % 局部峰值法
    p_fft = FFT(X1_eff);         % FFT方法
    
    if p > 0 && abs(p - p_fft) <= 1 && cv < 0.15
         % 若两种方法得到的周期较为接近且峰间隔波动较小，认为周期稳定
         switch round(p)
             case 1
                 catLabel = 11;  % P1
             case 2
                 catLabel = 10;  % P2
             case 3
                 catLabel = 9;   % P3
             case 4
                 catLabel = 8;   % P4
             case 5
                 catLabel = 7;   % P5
             case 6
                 catLabel = 6;   % P6
             case 7
                 catLabel = 5;   % P7
             otherwise
                 catLabel = 4;   % MP: 多周期
         end
    else
         % 若未检测到稳定周期，根据 E_sim 进一步划分
         if E_sim < 0.25
              catLabel = 3;  % QP: 准周期
         elseif E_sim < 0.35
              catLabel = 4;  % MP: 多周期
         else
              catLabel = 1;  % HC: 超混沌
         end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [p, cv] = Period(X)
% Period 利用局部峰值检测法判断时间序列的平均周期和波动程度
%
%   输出:
%       p  - 平均周期（取整），若无法检测则返回 0
%       cv - 峰间隔的变异系数 (CV = std/mean)
    [~, locs] = findpeaks(X);
    if length(locs) < 2
         p = 0;
         cv = Inf;
         return;
    end
    diffs = diff(locs);
    p = round(mean(diffs));
    cv = std(diffs) / mean(diffs);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p = FFT(X)
% FFT 利用快速傅里叶变换检测时间序列的主周期
%
%   输出:
%       p - FFT 检测得到的周期（以采样点为单位），若无法检测则返回 0
    N = length(X);
    Y = fft(X - mean(X));
    P2 = abs(Y/N);
    P1 = P2(1:floor(N/2)+1);
    P1(2:end-1) = 2*P1(2:end-1);
    [~, idx] = max(P1(2:end));  % 排除直流分量
    dominantFreq = (idx+1)/N;   % 主频率
    if dominantFreq == 0
         p = 0;
    else
         p = round(1/dominantFreq);
    end
end